Population structure and phylogeography of three closely related tree peonies

Abstract Paeonia decomposita, Paeonia rotundiloba, and Paeonia rockii are three closely related species of Sect. Moutan is distributed in the montane area of the Eastern Hengduan Mountain region. Understanding the population history of these three tree peony species could contribute to unraveling the evolutionary patterns of undergrowth species in this hotspot area. We used one nuclear DNA marker (internal transcribed spacer region, ITS) and two chloroplast DNA markers (matK, ycf1) to reconstruct the phylogeographic pattern of the populations. In total, 228 individuals from 17 populations of the three species were analyzed in this study. Three nuclear clades (Clade I – Clade III) and four maternal clades (Clade A – Clade D) were reconstructed. Molecular dating suggested that young lineages diverged during the late Pliocene and early Pleistocene, younger than the uplift of the Hengduan Mountains but older than the last glacial maximum (LGM). Significant population and phylogeographic structures were detected at both markers. Furthermore, the populations of these tree peonies were overall at equilibrium during the climatic oscillations of the Pleistocene. The simulated palaeoranges of the three species during the LGM period mostly overlapped, which could have led to cross‐breeding events. We propose an evolutionary scenario in which mountain orogenesis around the Hengduan Mountain area triggered parapatric isolation between maternal lineages of tree peonies. Subsequent climatic fluctuations drove migration and range recontact of these populations along the valleys. This detailed evolutionary history provides new insights into the phylogeographic pattern of species from mountain‐valley systems.


| INTRODUC TI ON
The East Asian flora is one of the most diverse temperate flora in the world (López-Pujol et al., 2011;Lu et al., 2018;Qi et al., 2012;Qiu et al., 2011;Wang et al., 2015). Flora from montane areas contributes to more than half of the plant diversity in East Asia, especially the mountain systems in middle to western China, such as the Qinling and Hengduan Mountains. In situ, diversification may have played a vital role in shaping the plant diversity of Hengduan forests, with geological history and climatic oscillations influencing the current distributions and genetic structures of many plant species (Favre et al., 2015;Perrigo et al., 2020;Xiang et al., 2017;Yang et al., 2022;Zhang et al., 2018). The Hengduan Mountains were thought to have experienced a series of topographical deformations from the late Oligocene (Ding et al., 2020;Su et al., 2019). Orogenic activities in mountain ranges increase topographic heterogeneity, create dispersal barriers or corridors, form novel habitat types, promote speciation, and influence phylogeographic patterns of species (Liu, Moller, et al., 2013;Yin et al., 2021). In addition, climate oscillations usually induce range contractions or expansions from refugial areas of many plants, with populations moving latitudinally along valleys (Tank & Sang, 2001;Vitt et al., 2010) and altitudinally on mountain slopes Zhang et al., 2018). The glacial-interglacial cycles of the Quaternary repeatedly altered organismal distributions and triggered unstable population dynamics (Matias et al., 2017).
Lineages of plant species may have opportunities to mix or to split due to isolation.
Lineage fusion and splitting were probably caused by orogeny and climate change. Lineage fusion is a special case in which two or more populations of parental lineages that were not reproductively isolated collapse into a single lineage after they meet, which has an important impact on present-day distributions of genetic diversity by reshuffling variation (Garrick et al., 2019(Garrick et al., , 2020Kearns et al., 2018).
Lineage fusion events can be realized by horizontal gene transfer or hybridization (Niehus et al., 2015;Yue, Hu, et al., 2012). Mountainvalley systems, a typical landscape of the Qinling-Hengduan Mountain region, are usually characterized by high mountains and deep valleys, and diverse local environments (Ding et al., 2020;Wang et al., 2022). If the orientation and morphology of mountain valleys are suitable, valleys and river drainages could serve as corridors for organisms around this area during glacial-interglacial cycles, promoting the occurrence of hybridization events and shaping the population structure of the plants here Du et al., 2022;Yue, Chen, et al., 2012;Yue & Sun, 2014;Zhang et al., 2011). In contrast, lineage splitting could be the most important mode of speciation when a mountain-valley acts as a natural geographic barrier (Yang et al., 2014;. However, researchers need more information on detailed lineage fusion and lineage splitting in mountain-valley ecosystems .  (Hong, 2021;Zhou, 2006). P. decomposita is mainly distributed along the Dadu River, and westernmost to Kangding, Sichuan; P. rotundiloba is mainly distributed along the upper reaches of the Minjiang River, and northernmost to Diebu, Gansu; and P. rockii is found in a relatively large area from central to western China (Thompson et al., 1994;Zhou et al., 2014). Previous studies showed that the three closely related tree peonies were morphologically different from each other (Hong, 2011b). P. decomposita is much more similar to P. rotundiloba, except for the differences in carpel number, disk coating degree, and leaflet morphology (Hong, 2011a). However, the phylogenetic relationships of the three closely related species are still controversial, as different genetic markers have emerged (Tank & Sang, 2001;Zhao et al., 2008;Zhou et al., 2014). Such discordance may indicate the effects of incomplete lineage sorting or introgressive hybridization (Vargas et al., 2017;Xu et al., 2019). Most populations of P. decomposita, P. rotundiloba, and P. rockii have been facing the risk of population extinction due to human disturbance, modern habitat fragmentation, and obstacles to seed reproduction (Song et al., 2012;Yuan et al., 2011). Therefore, more studies covering large numbers of populations on lineage history and population dynamics are needed to understand the population structure and speciation history and to formulate sustainable conservation strategies for these three species.
In this study, we present an analysis of the population genetic structure, phylogeography, and historical demography of P. rockii, P. decomposita, and P. rotundiloba using one nuclear DNA marker (internal transcribed spacer region, ITS) and two chloroplast DNA markers (matK, ycf1) and ecological niche modeling (ENM). Our aims in this study are to (i) reveal the geographic patterns of population genetic variation of the species; (ii) elucidate the evolutionary factors that might be responsible for the patterns and levels of genetic variation observed; and (iii) describe the demographic patterns of three closely related tree peonies that evolved during the LGM.

| Phylogenetic analysis
Sequences were initially subjected to a BLAST search against the GenBank database (www.ncbi.nlm.nih.gov) to detect contamination and confirm the targeted markers. Sequences were aligned using MAFFT v.7 and then adjusted manually in BioEdit v.7.0.0 (Hall, 1999).
All characters were weighted and unordered equally, and gaps were treated as missing data. The best nucleotide substitution models for cpDNA and nrITS were GTR and GTR + G, selected using the Akaike information criterion (AIC) in MrModeltest 2.3 (Nylander, 2004). For Bayesian analyses, two independent runs were performed through 10,000,000 generations with four Markov chains.

| Phylogeographic and population genetic data analyses
For both cpDNA (chlorotypes) and ITS, indices of haplotype diversity (h) and nucleotide diversity (π) were estimated according to Nei (1987) at the levels of populations (h S and π S ) and species (h T and π T ) using DnaSP (Librado & Rozas, 2009). Genealogical relationships among cpDNA and nDNA haplotypes identified were constructed using the program NETWORK v.5 (available at http:// www.fluxu s-engin eering.com/index.htm) with the MP criterion. A split phylogenetic tree based on the Neighbor-Net algorithm (NNet; Huson & Bryant, 2006) was also reconstructed for the nrITS haplotypes using SPLITSTREE v.4.14.6. Neighbor-Net networks can indicate any detailed potential conflicts among different haplotypes caused by complex evolutionary events such as hybridization, polyploidization, and recombination (Kilian et al., 2007;Liu et al., 2014). SAMOVA 1.0 (Dupanloup et al., 2002) was used to define the best groups of populations (K) that were geographically homogeneous and maximally differentiated from each other, with the number of groups ranging from 2 to 20. The analysis was reported five times for each K value with 1000 independent iterations, starting from 100 random initial conditions. The configuration with the highest and most significant F CT value was retained as the best grouping of populations.
Genetic diversity (H S ) and total genetic diversity (H T ) within a population, as well as genetic differentiation between populations (G ST and N ST ), were estimated using PERMUT v1.2.1 with the 1000 permutations test (http://www.pierr oton.inra.fr/genet ics/labo/ Softw are/Permu tCpSSR). The G ST and N ST were compared to determine the phylogeographic structure following Pons and Petit (1996).
Hierarchical analysis of molecular variance (AMOVA ;Excoffier et al., 1992) was performed to quantify the genetic differentiation among groups, within populations, and among populations using Arelquin v.3.0 (Excoffier et al., 2007), with 1000 random permutations. Populations were partitioned by geography or species.
Geographical groups were obtained from SAMOVA Correlations between the genetic distances (F ST ) and geographical distances of all populations (Mantel test) were tested using the software Arlequin v.3.0 (Excoffier et al., 2007) with 10,000 permutations at 0.05 significance levels. Pairwise geographical distance matrices were constructed by GenAlEx6 (Peakall & Smouse, 2012). Tajima's D (Tajima, 1989) and Fu's F S (Fu, 1997) were calculated in ARLEQUIN to test for deviations from neutrality. In an attempt to further infer demographic processes, we also tested the null hypothesis of spatial expansion and pure demographic expansion using mismatch distribution analysis (MDA) in ARLEQUIN (Excoffier et al., 2007) at the lineage and species levels. Goodnessof-fit was tested with the sum of squared deviations (SSD) between observed and expected mismatch distributions and Harpending's (Relethford & Harpending, 1994) raggedness index (HRag) using 1000 parametric bootstrap replicates for each model. Where the sudden expansion model was not rejected, the formula T = τ/2u (Rogers & Harpending, 1992) was used to estimate the age of expansion (t), where u = μkg; μ is the neutral mutation rate of the sequence in substitutions per site per year (s/s/y), k is the average sequence length (here, ITS, 636 bp; cpDNA, 1315 bp) and g is the generation time in years (Wang, 2010). Considering that there is no fossil record of Paeonia, a constant cpDNA substitution rate range for most angiosperm species was adopted (from 1.0 × 10 −9 s/s/y to 3.0 × 10 −9 s/s/y; Wolfe et al., 1987). For ITS, we adopted a nucleotide substitution rate ranging from 3.94 × 10 −9 to 6.06 × 10 −9 (Richardson et al., 2001).

| Demographic analyses and lineage divergence time of three closely related tree peonies
To estimate the change in demographic growth over the history of the samples, we carried out two Bayesian skyline plot analyses (BSPs), implemented in BEAST v1.8.2 (http://beast.commu nity/). The analyses were run using an uncorrelated relaxed molecular clock under Bayesian skyline models for 10 8 iterations with burn-in periods of 25%. The substitution models were set to GTR + I + G for nrITS and HKY + I for cpDNA. Genealogies and model parameters were sampled every 10,000 iterations. Chain convergence was inspected in Tracer v.1.5 software, and the results with an effective sample size (ESS) for each parameter >200 were accepted. The BSP was visualized in the program Tracer v.1.5, which summarizes the demographic history over time.
BEAST v.1.8.2 was also used to estimate divergence times since the most recent common ancestor of all the lineages. The program BEAUti (Drummond et al., 2012) was used to set criteria for the analysis, in which we applied the GTR + I + G substitution model for nrITS and the GTR + I substitution model for cpDNA, according to MrModeltest 2.3 (Nylander, 2004). Other nondefault settings were as follows: an uncorrelated relaxed molecular clock model (Drummond et al., 2002) was selected; the Yule Process was specified as the tree prior; and the prior parameters of "ucld.mean" were set to "uniform" distributions. Posterior distributions of parameters were obtained using MCMC analysis for 10 7 generations with 25% burn-in. The log file was then checked for convergence of the chains using Tracer 1.5.
A maximum clade credibility (MCC) tree with the maximum sum of posterior probabilities on its internal nodes  was adopted to summarize the samples from posterior distributions using TreeAnnotator v.1.8.2 , in which the posterior probability limit was set to 0.5 summarizing mean node heights. The MCC tree was visualized and edited using FigTree v.1.4.3, and the means and 95% higher posterior densities (HPDs) of the tree were obtained from it. The 95% HPD represents the shortest interval that contains 95% of the sampled values from the posterior (Drummond et al., 2006).

| Ecological niche modeling
To examine niche divergence between the targeted species, we used the maximum entropy model and machine-learning algorithm implemented in MAXENT v.3.3.3k Phillips & Dudík, 2008) to predict ecological niche models based on presenceonly data obtained from our field observations, herbarium records in the Chinese Virtual Herbarium (http://www.cvh.ac.cn/) and the literature (Yuan et al., 2012) for all three species (Table S6). Current values for 19 bioclimatic variables together with the altitude data package, which was at a resolution of 2.5 (arc-minutes), were downloaded from the WorldClim dataset set (Hijmans et al., 2005) as environmental layers. All 228 individuals were used for ENM analysis, and the number of individuals of each species was above 50. We

| RE SULTS
In this study, 221 individuals were successfully sequenced for cpDNA, and 228 individuals were sequenced for nrITS (Tables S1 and S3). The total length of concatenated chloroplast DNA alignment was 1315 bp, containing 8 nucleotide substitutions and 4 indels (insertion/deletion).

| Phylogenetic analyses and estimations of divergence times
Four maternal clades were reconstructed (Clade A-D, Figure 4) based on chlorotype data, in which the paraphyly of P. decomposita and P. rotundiloba was confirmed in this analysis using the neighbor-joining (NJ), Bayesian, and maximum likelihood methods. Clade C contained haplotypes from P. decomposita and P. rotundiloba. P. rockii was polyphyletic and included haplotypes C7-C9 (extracted from populations JS, WD, WQ) and C10-C13 (extracted from population CZ), corresponding to Clades B and D respectively (for details, see Table S1 and Figure 4). Four clades were connected by one median vector in the phylogenetic network (Figure 1c). For the nrITS database, the NJ tree, MrBayes tree, ML tree, and Splitstree indicated minor differences.
Three main clades were recovered (Clades I-III). While Clades I and III harbored haplotypes extracted from both P. decomposita and P. rockii, Clade II contained haplotypes from all three tree peonies (for details, see Table 2 the cpDNA chronogram, BEAST adopted substitution rates ranging from 0.3 × 10 −9 to 1.0 × 10 −9 substitutions per site per year (s/s/y) (Su et al., 2015). For the nrITS chronogram, the substitution rate was defined from 3.94 × 10 −9 to 6.06 × 10 −9 substitutions per site per year (s/s/y) (Richardson et al., 2001). We estimated the crown age of three lineages and outgroups as c. 2.5651 Ma (95% HPD: 1.3941-3.929 Ma; node 1), which fell into the late Pliocene to the early Pleistocene. The split of lineage III was estimated to be c. 2.2029 Ma (95% HPD: 1.2214-3.2865 Ma; node 2), and the lineage split between lineages I and II was dated to c. 2.0463 Ma (95% HPD: 1.1751-3.0404 Ma; node 3).  Table S1); the solid black lines represent three groups according to SAMOVA, and the dotted lines display six groups (b). Maximum parsimony network based on chlorotypes (c). The size of the pie is positively correlated with the frequency of each haplotype, and black dots represent haplotypes that are missing, extinct, or not observed. haplotypes (see Table S2 for details), which could be clustered into three lineages (red 49.64%, blue 10.95%, and purple 39.42%). The

| Phylogeographic and population structures
Central groups of the two markers were both located mainly along the valleys of the Min-Jiang River in northwestern Sichuan Province. AMOVA revealed that the differences among populations (for total populations) explained 91.7% of the total cpDNA variation.
However, variation within populations accounted for most (83.68%) of the variation in nrITS (Figures S2 and S3). For the three groups identified by SAMOVA, AMOVA showed that a large proportion (52.86%) of the cpDNA variation was partitioned among groups, suggesting a significant regional population substructure. In contrast, only a small proportion (16.11%) of the nrITS variation was distributed among groups (Tables 1 and 2).

| Demographic analyses
The null hypotheses of population expansion were rejected. No sig-  and showed a slight decrease during the late Pleistocene ( Figure S4), while the effective population sizes based on nrITS were shown to be slightly increased during the middle Pleistocene, albeit along with a significant increase approximately 120,000 years ago ( Figure S5).

| Present and past (LGM) distribution
The AUC values for the current potential distributions of P. decomposita, P. rotundiloba, and P. rockii indicated good predictive model performance (0.990, 0.995, and 0.992, respectively). The current distributional predictions covered the species' extant distribution.
However, there were also some predicted areas where the species are not found at present. In contrast, distributional predictions for the LGM (c. 21,000 year BP) were different from the present, revealing a general southwards range shift. Furthermore, the simulated distribution areas of the three tree peonies at the LGM partially overlapped ( Figure 6).

| Genetic diversity and phylogeographic pattern
The three tree peonies are characterized by high genetic diversity at the overall population level in both cpDNA and nrITS markers. In this study, the h T values of cpDNA and nrITS were as high as 0.835 and 0.943, respectively, and were significantly higher than those F I G U R E 3 Maximum parsimony network based on nrITS haplotypes. The size of the pie is positively correlated with the frequency of each haplotype, and black dots represent haplotypes that are missing, extinct, or not observed. Lineage Ⅰ Lineage Ⅲ Lineage Ⅱ Outgroup reviewed by Qiu et al. (2011). The three tree peonies mainly propagate by seeds (Cheng et al., 1997), which are normally dispersed by either gravity or rats (Hong, 2011b). Individuals always start to flower 4-6 years after germination, and it is generally assumed that the generation time of Paeonia species is 10 years (Sang et al., 1997;Xu et al., 2019). The long generation time of tree peonies may have acted as an intrinsic buffer against the loss of genetic diversity because individuals of different generations could share the same habitat (Hailer et al., 2006). Individual tree peonies could live for more than 30 years according to our long-term field observations by F I G U R E 4 BEAST-derived chronograms based on cpDNA haplotypes. BEAST-derived topology is similar to Bayesian/ML/MP analyses. C1-C13 represent the ingroups; details are referenced in Table 1. C14-C20 represent the outgroups; C14, P. lutea; C15-C19, P. devalayi; C20, P. ludlowii.   which could accumulate haplotype diversity in the maternal system (Chen et al., 2017;Montarry et al., 2019). For nrITS markers, tree peonies are pollinated by insects such as bees, and pollen flow via insects can make tree peonies exhibit a more uniform distribution of diversity even in sparse Cupressus chengiana forests, young secondary deciduous broad-leaved forests (Wang, 2020;Xu et al., 2019;Yuan et al., 2011). In addition, hybridization, either cross-species or intraspecific hybridization, is increasingly recognized as a generator of diversity (Lamichhaney et al., 2018;Li et al., 2020). Our results were consistent with several previous studies in the Hengduan Mountains, such as those studying Ligularia tongolensis (Wang et al., 2011), Terminalia franchetii , Cyananthus delavayi (Li et al., 2012), Lespedeza buergeri (Jin et al., 2016), and Forsythia suspense (Li et al., 2011).
The three species of tree peonies showed significant genetic differentiation and phylogeographic structure within distinct topo- serve as a geographical barrier, as steep mountains and thick forest vegetation prevent seed and pollen dispersal from crossing, to realize gene transfer between tree peony recipients and potential donors.
Different levels of genetic differentiation and genetic structure among populations between cpDNA and nrITS could be attributed to the characteristics of their different inheritance systems. The cpDNA is maternally inherited with low recombination and mutation rates and is free of the types of phenomena that lead to violation of crossbreeding (Mark et al., 2007;Schaal et al., 1998). In contrast, the nuclear rDNA ITS is parentally inherited with relatively higher recombination and substitution rates (Mark et al., 2007). Moreover, insect-triggered pollen flows carrying genes, such as ITSs, are always stronger and farther from the original individuals than the seed flows, which may increase the gene introgression between the populations, lead to the unification processes of alleles and form a lower phylogeographic structure across the range.

| Lineage sorting and fusion of the three Paeonia species
Complete lineage sorting was observed from the chloroplast segments. (Figures 4 and 1c). The splits of the main lineages were dated from the late Pliocene to the early Pleistocene (c. 3.23 Ma to c. 2.15 Ma) for the cpDNA markers (Figure 4), both of which were TA B L E 1 The analysis of molecular variance (AMOVA) for cpDNA data. forming a mountain-valley landscape in this and adjacent areas (Clark et al., 2005;Wang et al., 2012;Xing & Ree, 2017) and inducing isolation and allopatric speciation events. Close correlations between lineage splitting and geologic changes have also been observed in other plants, e.g., Tetrastigma hemsleyanum (Wang et al., 2015), Taxus wallichiana (Liu, Moller, et al., 2013), Picea purpurea ,

Sum of squares
Hippophae tibetana , and Picea likiangensis . In this study, the uplift of the Hengduan Mountains could also lead to the formation of new lineages of the three tree peonies.
From the nuclear genome data, lineage fusion was detected among the three tree peonies (Figures 5, 2c and 3). There are several reasons for lineage fusion (Yue, Hu, et al., 2012), and hybridization events could more reasonably explain the fusion of the lineages in this study. Due to the concerted evolution of intra-and interchromosomal loci in the long-term process of evolution, multiple copies of the nuclear ribosomal ITS (nrITS) within a genome tend to highly homogenize the variation among nrITS repeats, as the nrITS can be found inside individuals of one species (Vollmer & Palumbi, 2004).
Lineages of nrITS were split prior to the early Pleistocene ( Figure 5), which coincides with the intense uplift of the QTP. From the early Pleistocene to the present, the rate of concerted evolution is sufficient to make the ITS remain a single homogeneous copy (Baldwin et al., 1995;Leitch et al., 2004). Older lineages occurring in younger individuals and populations indicate hybridization events, with perhaps more than one event, leading to the lineage fusion pattern of the three tree peonies.

| Population dynamics during the Quaternary leading to lineage fusion
Lineage fusion could be realized by population expansions or migration induced by climate and geologic changes. Demographic or range expansions may lead to short-distance pollen flow exchange, as well as massive gene introgression of different species (Excoffier et al., 2009;Niehus et al., 2015). We performed demographic and ecological niche modeling analyses to determine lineage fusion.
However, multiple lines of evidence from the neutrality test, mismatch distribution analysis, and BSP analysis presented here suggested that range expansions or contractions did not occur during the LGM (Figures S4 and S5).
However, the Quaternary climatic fluctuation and plateau uplift have largely shaped the current patterns of species diversity (Xing & Ree, 2017;Yu et al., 2019). Climate oscillations impel latitudinal/altitudinal shifts in species distribution ranges, which may increase the chance of hybridization between species at their boundaries (Garroway et al., 2010;Mimura et al., 2014). Because of the extremely cold climate during the LGM, the three closely related peonies migrated to the lower plate through the valleys between mountains ( Figure 6). Therefore, the ranges of these three species probably overlapped at least in part (in the Chengdu Plain according to the current geographical division), which would have provided an opportunity for short-distance pollen dispersal and further hybridization (Schweiger et al., 2008).
Insect pollinators are always extremely efficient. Pollen flows carried by insects could unify the ITS copies among the three species as soon as possible once the metapopulation around the Chengdu Plain area occurred. We performed demographic and ecological niche modeling analyses to determine lineage fusion. Ecological niche modeling showed that the simulated core distribution areas during the LGM were not far from the present distribution areas for each tree peony ( Figure 6). Similar results have been observed in Europe, North America, and the middle and eastern Tibetan Plateau, revealing that long-distance migration from old core distribution areas to new ones can follow climatic oscillations (Clayton et al., 2009;Givnish et al., 2016;Liu, Ickert-Bond, et al., 2013;Nakamura et al., 2012).

| Conservation strategies for three closely related species
Due to the frequent occurrence of geological disasters in the distribution of peony in Southwest China, many wild resources have been destroyed. The pollination period coincides with the rainy season, which reduces insect pollination efficiency. Moreover, with the promotion of peony oil in recent years, the picking of wild germplasm resources has become more frequent, which is not conducive to the natural renewal of wild populations and results in smaller population size. Small-scale populations do not exchange genes with the outside world and are more likely to face the risk of genetic drift. In the long run, the plant resources of the peony group are bound to face the danger of extinction. In addition, this study shows that the populations of the three related species of the tree peony group contain abundant haplotypes, and the haplotypes contained in each species are highly endemic. Therefore, we believe that existing populations can be protected in situ as protection units, and these populations can be classified into nature reserves to reduce the effects of human disturbance. Second, we encourage research on the mechanism of peony endangerment. Finally, relocation protection can be carried out. Our field investigations found that peony regeneration is better in cases of less human disturbance, which can provide better conditions for peony reproduction through introduction and conservation.

| CON CLUS IONS
Our study of P. rockii, P. decomposita, and P. rotundiloba, integrating molecular phylogeography, and ecological niche modeling, strongly suggests that the observed patterns of genetic variation and evolutionary history are best explained by a combination of late Miocene and late Pliocene geological activity and late Quaternary climatic fluctuation.
Furthermore, valleys among mountains acted as corridors, playing important roles in lineage confluence and species migration.
The glacial periods could have triggered the populations to migrate The results of this study will increase our comprehensive understanding of the origin, evolution, phylogenetic relationship, and distribution of sect. Moutan. According to our results, we propose that these three tree peonies are not fully fledged sister species but instead are geographical races, semispecies, or subspecies that are undergoing an allopatric phase of divergence. However, due to the evaluation of limited samples and only a few molecular markers, these results may need validation through additional research.
In a future study, multilocus approaches using many genes, or even a whole-genome sequencing approach, should be applied to resolve the species status of these three closely related tree peonies.

ACK N OWLED G M ENTS
We thank Hao Wang, Yong Yang, Tao Zhong, and Yinyin Jiang for their kind help in field assistance. Special thanks go to Kangshan Mao (Sichuan University) for his constructive comments on experiment design and Ecological niche modeling analysis. We would like to thank Stephen Gaughran at Yale University for his assistance with the English language and grammatical editing of the manuscript.

FU N D I N G I N FO R M ATI O N
None.

CO N FLI C T O F I NTER E S T S TATEM ENT
The authors declare no conflicts of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
The ITS sequences of per individual of have been deposited in GenBank P. rotundiloba and P. decomposita (accession numbers OM570335-OM570549), P. rockii (accession numbers OM562281-OM562376).
The ycf1 sequences of per individual have been deposited in GenBank (accession numbers Om993806-om994027).
The matK sequences of per individual have been deposited in GenBank (accession numbers Om993584-om993805).